CONVERGENT FINITE DIFFERENCE METHODS FOR 
ONE-DIMENSIONAL FULLY NONLINEAR SECOND ORDER 
PARTIAL DIFFERENTIAL EQUATIONS 

XIAOBING FENG*, CHIU-YEN KAQt, AND THOMAS LEWIS* 

Abstract. This paper develops a new framework for designing and analyzing convergent finite 
difference methods for approximating both classical and viscosity solutions of second order fully 
nonlinear partial differential equations (PDEs) in 1-D. The goal of the paper is to extend the successful 
framework of monotone, consistent, and stable finite difference methods for first order fully nonlinear 
Hamilton-Jacobi equations to second order fully nonlinear PDEs such as Monge-Ampere and Bellman 
type equations. New concepts of consistency, generalized monotonicity, and stability are introduced; 
among them, the generalized monotonicity and consistency, which are easier to verify in practice, 
are natural extensions of the corresponding notions of finite difference methods for first order fully 
nonlinear Hamilton-Jacobi equations. The main component of the proposed framework is the concept 
of a "numerical operator", and the main idea used to design consistent, generalized monotone and 
stable finite difference methods is the concept of a "numerical moment" . These two new concepts 
play the same roles the "numerical Hamiltonian" and the "numerical viscosity" play in the finite 
difference framework for first order fully nonlinear Hamilton-Jacobi equations. In the paper, two 
classes of consistent and monotone finite difference methods are proposed for second order fully 
nonlinear PDEs. The first class contains Lax-Friedrichs-like methods which also are proved to be 
stable, and the second class contains Godunov-like methods. Numerical results are also presented 
to gauge the performance of the proposed finite difference methods and to validate the theoretical 
results of the paper. 
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1. Introduction. Fully nonlinear partial differential equations (PDEs) refers to 
a class of nonlinear PDEs which are nonlinear in the highest order derivatives of the 
unknown functions appearing in the equations. For example, the general first and 
second order fully nonlinear PDEs, respectively, have the form H{Wu,u,x) = and 
F{D^u, Vu, u, x) — 0, where Vu and D^u denote the gradient vector and Hessian 
matrix of the unknown function u. Fully nonlinear PDEs, which have experienced 
extensive analytical developments in the past thirty years (cf. [H [U [TT] [IH]), arise 
from many scientific and engineering applications such as differential geometry, astro- 
physics, antenna design, image processing, optimal control, optimal mass transport, 
and geostrophical fluid dynamics. Fully nonlinear PDEs play a critical role for the so- 
lutions of these applications because they appear one way or another in the governing 
equations of these problems. 

As expected, the study of first order fully nonlinear PDEs came first. Since the 
introduction of the notion of viscosity solutions by Crandall and Lions [6 in 1983, 
the past thirty years has been a period of explosive developments in analyzing first 
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order fully nonlinear PDEs. Starting with the pioneering work of Crandall and Lions 
[7j, extensive research has also been successfully carried out on developing numerical 
methods, in particular monotone as well as other types of finite difference methods, 
for computing viscosity solutions of first order fully nonlinear PDEs, especially those 
arising from the level set formulations of moving interfaces and those arising from opti- 
mal control (cf. [20] and the references therein). To overcome the low order accuracy 
barrier of monotone finite difference methods, various high order local discontinu- 
ous Galerkin (LDG) methods have also been developed recently in the literature (cf. 
PPI and the references therein) . 

In contrast with the success of PDE analysis and numerical approximation for 
first order fully nonlinear PDEs, the situation for second order fully nonlinear PDEs 
is very different. On one hand, like in the case of first order fully nonlinear PDEs, 
tremendous progresses in PDE analysis have been made in the past thirty years (cf. 
[mil]). On the other hand, not much progress on developing accurate and efficient 
numerical methods, especially Galerkin-type methods, for second order fully nonlinear 
PDEs has been made until very recently (cf. [13 [H] and the references therein) . The 
lack of progress is mainly due to the following two facts: (i) the notion of viscosity 
solutions is nonvariational; (ii) the conditional uniqueness (i.e., uniqueness only holds 
in a restrictive function class) of viscosity solutions is difficult to handle at the discrete 
level. The first difficulty prevents a direct construction of Galerkin-type methods and 
forces one to use indirect approaches as done in [HI dHl [HI [11] for approximating 
viscosity solutions. The second difficulty prevents any straightforward construction 
of finite difference methods because such a method does not have a mechanism to 
enforce the conditional uniqueness and often fails to capture the sought-after viscosity 
solution. Since the scope of this paper is confined to the finite difference method, 
Galerkin-type methods will not be discussed here. We refer the reader to the review 
paper [T3] for a detailed discussion of recent developments on Galerkin-type methods 
for second order fully nonlinear PDEs. 

The primary goal of this paper is to develop a new framework for designing and 
analyzing convergent finite difference methods for second order fully nonlinear (ellip- 
tic) PDEs. For the ease of presenting the ideas and to observe the page limitation of 
the journal, we shall only consider one-dimensional PDEs in this paper and leave the 
high dimensional generalizations to a forthcoming companion paper [TS]. We use the 
phrase "new framework" to distinguish the framework of this paper from the exist- 
ing (abstract) framework originally developed by Barles and Souganidis in |2j twenty 
years ago and further developed recently by Caffarelli and Souganidis in |5]. Unlike 
Barles and Souganidis' framework which is abstract and broader in applications, our 
framework is specifically and only designed for finite difference methods which can 
be easily implemented on computers. As a result, the proposed framework has the 
advantages of being simple to understand and easy to utilize in practice. Moreover, 
the new framework is a natural extension of the successful monotone finite difference 
framework developed for first order fully nonlinear Hamilton- Jacobi equations (cf. 
[3 [10] and the references therein) . The main concept of the new framework is the 
"numerical operator" . The key components of the framework are new and easy-to- 
check notions of consistency and generalized monotonicity (g-monotonicity) , which 
together with the well-known notion of stability, form the backbones of the proposed 
finite difference framework. After the framework is established, one must address a 
harder question of how to construct specific finite difference methods which fulfill the 
structure conditions (i.e., consistency, g-monotonicity, and stability) of the framework 
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in order to make the framework practically useful. Wc note that this question was 
not addressed in [2 as the goal of that paper was not to develop practical numerical 
methods, and it took seventeen years to construct the first finite difference method 
which fulfills the structure conditions laid out in [2] for the second order fully nonlin- 
ear Monge- Ampere equation in [19]. Moreover, the method of [19] is a nonstandard 
finite difference method because it requires the use of wide-stencil grids. We do want 
to remark that many numerical methods, which may or may not fulfill the structure 
conditions of [2], have been developed for Bellman type equations (cf. [31 Uni [13] 
and the references therein). To address the above key question, our main idea is 
to introduce a new concept called the "numerical moment". We like to stress that 
the numerical moment not only helps the construction of desired g-monotone finite 
difference methods, but also, we believe, provides a fundamental and indispensable 
mechanism for a finite difference method to overcome the two major difficulties asso- 
ciated with numerical approximations of second order fully nonlinear PDEs. We also 
note that the new concepts of "numerical operators" and "numerical moments" for 
second order fully nonlinear PDEs are natural extensions of the well-known concepts 
of "numerical Hamiltonians" and "numerical viscosities" for first order fully nonlinear 
Hamilton- Jacobi equations. 

This paper is organized as follows. In Section [2] we collect some preliminary 
materials such as notation and definitions. In Section [3] we present our finite dif- 
ference framework. The motivation and main ideas are heuristically explained. The 
main concepts and definitions of numerical operators, consistency, g-monotonicity, 
and stability are formally introduced and defined. The main result of this section is 
a convergence theorem which asserts that the solution of any consistent, g-monotone 
and stable finite difference method is guaranteed to converge to the unique viscosity 
solution of the underlying second order fully nonlinear PDE. In Section |4| we intro- 
duce the concept of a numerical moment. With the help of the numerical moment 
and the inspiration given by the convergent finite difference schemes for first order 
fully nonlinear Hamilton-Jacobi equations, we are able to construct two classes of 
consistent and g-monotone finite difference methods. The first class contains Lax- 
Friedrichs-like methods and the second class contains Godunov-like methods. By 
using a non-standard fixed point argument we also prove that every consistent and g- 
monotone Lax-Friedrichs-like method is uniquely solvable and stable for a given class 
of fully nonlinear operators. In Section [5] we present some detailed numerical results 
to gauge the performance of the proposed finite difference methods and to validate 
the theoretical results of the paper. The paper is concluded by a short summary in 
Section |6l 

2. Preliminaries. In this paper we adopt standard function and space notations 
as in [TTJ [3]. For example, for a bounded open domain H. C R'', B{n), USC{fl) 
and LSC{i}) are used to denote, respectively, the spaces of bounded, upper semi- 
continuous and lower semicontinuous functions on fl. Also, for any v € -6(0), we 
define 

V* (x) := limsupw(j/) and v^^,{x) := liminf ?;(?;). 

y — y 

Then, v* € USC{Q) and t;* € LSC{Q), and they are called the upper and lower 
semicontinuous envelopes of v, respectively. 

Given a bounded function F : S'^^''' x R"^ x R x 17 R, where S"^^^ denotes the 
set of d X d symmetric real matrices, the general second order fully nonlinear PDE 



4 



X.FENG, C.-Y.KAO, AND T. LEWIS 



takes the form 

(2.1) F{D'^u,Wu,u,x) = inU. 

Note that here we have used the convention of writing the boundary condition as a 
discontinuity of the PDE (cf. ^ p.274]). 

The following two definitions are standard (cf. pTl HI 1^). 

Definition 2.1. Equation ( |2.1| ) is said to be elliptic if for all (p, A,a;) € R'' x 
R X r2 there holds 

(2.2) F{A,p,X,x) < F{B,p,X,x) yA,BeS'^'"^,A>B, 



where A> B means that A — B is a nonnegative definite matrix. We note that when 

OF 
dA 



F is differentiable, the ellipticity also can be defined by requiring that the matrix 



is negative semi-definite (cf. [TTJ p. 441]). 

Definition 2.2. A function u G B(yi) is called a viscosity suhsolution (resp. 



supersolution) of (2.1 1 if, for all Lp S C^(f2), if u* — ip (resp. u^, — Lp) has a local 



maximum (resp. minimum) at xq £ i^, then we have 

F^{D'^ip{xo),W(p{xo),u*{xo),Xo) <0 
F*{D'^lp{xq),'V(p{xo),u^{xo),xq) > 0). The function u is said to be a vis- 



cosity solution of (2.11 if it is simultaneously a viscosity suhsolution and a viscosity 
supersolution of ( |2.1[ ). 

We remark that if F and u are continuous, then the upper and lower * indices 
can be removed in Definition |2.2[ The definition of the ellipticity implies that the 
differential operator F must be non-increasing in its first argument in order to be 
elliptic. It turns out that the ellipticity provides a sufficient condition for equation 



(2.1) to fulfill a maximum principle (cf. [TTl|4|). It is clear from the above definition 
that viscosity solutions in general do not satisfy the underlying PDEs in a tangible 
sense, and the concept of viscosity solutions is nonvariational. Such a solution is 
not defined through integration by parts against arbitrary test functions; hence, it 
does not satisfy an integral identity. As pointed out in Section [T] the nonvariational 
nature of viscosity solutions is the main obstacle that prevents direct construction of 
Galerkin-type methods, which are based on variational formulations. 

3. A monotone finite diflference framework. We consider the following fully 
nonlinear second order two-point boundary value problem: 

(3.1) F{uxx,x)=Q, a<x<b, 

(3.2) u{a) = Ua, 

(3.3) u{b) = Ub, 

where Ua and Ub are two given numbers and F is assumed to be an elliptic operator 
in a function class A C C°(S1). We remark that the results of this paper can be easily 
extended to PDEs with general form F{u,j,x,u,j,,u,x) = 0. 

To construct finite difference methods for the above problem, we first need to 
have a mesh for the domain/interval f2 := {a,b). For simplicity, we only consider 
uniform meshes here, although our methods can be easily generalized to nonuniform 
meshes. Let J be a positive integer and h = jEf • We divide ft into J — 1 subin- 
tervals/subdomains with grid points Xj = a + (j — l)h for j = 1, 2, • • • , J, and let 
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Th — {xjYj^i be a mesh of ^l. Define the forward and backward difference operators 

by 

, v(x + h)-v{x) v(x) - v(x — h) 
djv{x) , 5^ v{x) := , 



for a continuous function v defined in fl and 



for a grid function V defined on the mesh T/i- The operators and b~ will serve as 
building blocks in the construction of our finite difference methods in the sense that we 
approximate all first and second derivatives by using combinations and compositions 
of these two operators. 

To approximate u^ixj), we have two options 

As a result, we have three possible ways to approximate ^^^.(xj) given by 
Uxx(x^) « btbtu{xj), u^xixj) « 8-8-u{x^), 

ft is easy to verify that 

(5+(5+u(a;j) = (5^u(a;j-+i), (5^(5^w(xj) = (5^w(xj_i), (5^(5+u(a;j) = 5lu{xj), 
where 

v(x ^ h) — 2v{x) + v{x + h) 



5^v{x) := 



for a continuous function v and 



for a grid function V on the mesh Th- 

The above simple argument motivates us to propose the following general finite 
difference method for equation (3.f ): Find a grid function U such that 

(3.4) F{6lU,^uSlU„SlU,+^,x,) = 

for j = 2, 3, • • • ,J — 1. As expected, Uj is intended to be an approximation of u{xj) 
for j = 1, 2, • • • , J, and Uq and t/j+i are two ghost values. 

Definition 3.f. The function i^:RxRxRxR-^Rm (3.4| is called a 
numerical operator. Finite difference method (3.4 1 is said to be an admissible scheme 
for problem (3.f )-(3.3l if it has at least one (grid function) solution U such that 
Ui — Ua and Uj — ujj. 

It is easy to understand that F needs to be some approximation of the differential 
operator F in order for scheme (3.4) to be relevant to the original PDE problem. 
Generally, different numerical operators F should result in different finite difference 
methods. A natural and important question is how to construct F . We shall defer 
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answering this question to the next section where we present two types of numerical 
operators F. For now, we propose a set of conditions (or properties) which we hke to 
impose on F. We choose conditions such that if F satisfies them, then the solution of 



the finite difference method (3.4 1 is guaranteed to converge to the viscosity solution 



of problem (3.1)-(3.3). The conditions will be reflected in the following definition. 



Definition 3.2. 

(i) Finite difference method (3.4) is said to be a consistent scheme if F satisfies 



(3.5) 


lim inf 

p;.— *-p,fc=l,2 


3 Fipi,P2,P3,0 > F*{p,x), 






(3.6) 


lim sup 


Fip,,p2,P3,0<F*{p,x), 

3 




— >-p 1 fc= 1 , 2 




(3.7) 


lim inf 

— >■ — oo , fc = 1 , 2 


^ F{pi,P2,P3,0 > ^'*(-oo,a;) 






(3.8) 


lim sup 

Pfc^oo,fc = l,2 
5-+X 


F{pi,P2,P3,£.) < F*{oo,x) : 

3 



for p G R, where F^, and F* denote respectively the lower and the upper 
semi- continuous envelopes of F. 



(ii) Finite difference method (3.4) is said to be a g-monotone scheme if for each 
2 < J < J — 1, F(pi,p2,P3,Xj) is monotone increasing in pi and ps and 
monotone decreasing in p2; that is, -F(t, J,, f, Xj) for j = 2,3, ••• , J — 1. 



(iii) Let (3.4) be an admissible finite difference method. A solution U of (3.4) is 
said to be stable if there exists a constant C > 0, which is independent of h, 
such that U satisfies 

(3.9) IMi^^n) := max^|C/,|<C. 



Also, ( |3.4[ ) is said to be a stable scheme if all of its solutions are stable 
solutions. 

Remark 1. (a) The consistency and g-monotonicity (generalized monotonicity) 
defined above are different from those given in \17[ F is asked to be monotone 
in S'^Uj-i,S'^Uj and S'^Uj+i, not in each individual entry Uj. To avoid confusion, we 
use the words "g-monotonicity" and "g-monotone" to indicate that the monotonicity 
is defined as above. We shall demonstrate in the next section that the above new 
definitions, especially the one for g-monotonicity, are more suitable and much easier 
to verify for (practical) finite difference methods. The new notions of consistency and 
g-monotonicity are logical extensions of their widely used counterparts for the first 
order Hamilton- Jacobi equations \2U^ . 

(b) On the other hand, the above stability definition is the same as that given in 

li[73[3'. 

(c) We note that if F is a continuous function, we can also assume that F is 
a continuous function. Then, (3.5) and (3.6) reduce to the condition F{p,p,p,x) = 



F{p,x). 

(d) The "good" numerical operators F we construct so far (cf. Section^ all have 
the form 

(3.10) Fip,,p2,P3,O^G{p^,P2,0 
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for some function G and p2 {pi +p^)/2. In other words, F is a function of p2 and 
P2 . Hence, a g-monotone F should be increasing in pi + p^ and decreasing in p2 ■ In 
this case, the consistency condition reduces to 



(3.11) 
(3.12) 
(3.13) 
(3.14) 



lim inf G{ai ,(T2,C) > F*{'P,x), 
lim sup G{ai ,a2,0 < F*iP>x)^ 

CTl ,cr2 -^P 
{-»x 

liminf G{(Ti,cr2,£,) > F^,{~oo,x) 



^1 •"2- 
£->x 



lim inf F(j),x), 



lim sup G(o'i, 0^2, < F* (oo, x) :— lim sup F{p, x) 



We shall need to use the above form of F in the proof of our convergence theorem, see 



Theorem 3.4 below. 



For a given grid function U, we define a piecewise constant extension function Uh 
of U as follows: 



(3.15) 



Uh{x) :— Uj Vx G 



] 



1,2, 



where 



±1 for j = l,2,... ,J. 



Definition 3.3. Problem (3.1|-(3.3) is said to satisfy a comparison principle if 



the following statement holds. For any upper semi-continuous function u and lower 
semi-continuous function v on fi, if u is a viscosity subsolution and v is a viscosity 



supersolution of (3.1)-(3.3|, then u<v onVt. 

Remark 2. Since the comparison principle immediately infers the uniqueness of 



viscosity solutions, it is also called a strong uniqueness property for problem (3.1) 



Q (cf W- 

We are now ready to state and prove the following convergence theorem, which 
is the main result of this paper. 

Theorem 3.4. Suppose problem (3.1)-(3.3) satisfies the comparison principle of 
Definition \3.3\ and has a unique continuous viscosity solution u. Let U be a solution to 
a consistent, g-monotone, and stable finite difference method (3.4 1 with F satisfying 
(3.10), and let u/j be its piecewise constant extension as defined above. Then u/i 



converges to u locally uniformly as h ^ 0^ . 
Proof. We divide the proof into five steps. 

Step 1: Since U satisfies (3.9), it is trivial to check that Uh satisfies 



(3.16) 

Define u,m G L°^{Vl) by 



\\uh\\L=°{n) < C. 



M(a;) := limsupu;i(^), u(a;) 

£->x 



lim inf Uh{S,)- 

h^0 + 



We now show that u and u are, respectively, a viscosity subsolution and a viscosity 



supersolution of (3.1)-(3.3). Hence, they must coincide by the comparison principle. 

Suppose that u—ip takes a local maximum at Xq G for some Lp e C^(i7). We first 
assume that </? e P2, the set of all quadratic polynomials. In Step 3 we will consider 
the general case ip G C^(ri). Without loss of generality, we assume u{xq) — <pi{x[)) 
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is a strict local maximum and u{xa) — ip{xq) (after a translation in the dependent 
variable). Then there exists a ball/interval, Br„{^a)i centered at xq with radius tq > 
such that 

(3.17) u{x) ~ Lp{x) < u{xq) ~ tf{xi^) = Q Vx e Broixo). 

Thus, there exists sequences {hk}k>i and {^k}k>i such that as fc — > oo, 
hk 0+, -> 2^0, ^i/ifc(Cfe) u{xo), 

u/jj.(x) — (p{x) takes a local maximum at ^k for sufhciently large k, 

and 
(3.18) 

where 



lim S^j^ Uh^i^k) =liminf(5^^u(a;o), 



UhiC - p) - 2uhiO + + P) 



P^ 



V^e {a + p,b~p), p>0. 



We remark that the right-hand side of (|3.18|) could either be finite or negative infinite. 
Tl 

(3.19) 



Then, there exists ko » 1 such that h}^ < and 

k^oo 



i — u/!.fc(&) - ^i^k) > Uhkix) - v{x) yx e Braixa), k > fco- 



Step 2: Since U satisfies (3.4 1 with F being of the form (3.10) at every interior 
grid point, it is easy to check that for x e il/j := (a + & — ^^JJ 

(3.20) = F{Sl ,^Uh{x - h), dl f^Uhix),Sl f^Uh{x + h),x) 

= G{SlM'^h{x),Sl,^Uk{x),x), 

SlMUhix) := Slj^Uhix -h) + Slj^Uhix + h). 
Since Uh^ix) — ^{x) takes a local maximum at and hk < ro for A; > fco, by 



where 



(3.19) we have 



(3.21) 



Sl,h^uh,{£,k) < Sl^hM^k) = <fixx{xo) Vfc > fco. 



Also, by ( |3.17| ), we get 

SImH^o) < SljM^o) = fxxixo) Vft. < ro. 

Thus, 

(3.22) limsup^^ hH^o) < ^fixxixo)- 

/i-i-O 

Next, a direct computation yields that 

(3.23) 5lj^Uh{x) ^ Sl f^Uh{x) + 2RhU,i{x), 
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where 



RhUh{x) := 5l 2hUh{x) - 5l uUh{x). 



By (3.181 and the definition of liminf we get 



(3.24) 



Thus, 



Hminf ,5^ 2/ifc"/iJf/c) = liminf (5^ 2hfe"(a;o) 

fe— >-oo fc-+oo 



/i— »-0 fc— >oo 



Uminf ^ liminf (5^ 2^,-^/1 J^fc) " ,1™ ^xM^hk^S.k) > 0, 

k^oo k^oo k-~^oo 



(3.25) 

and there exists a sequence {efc}fc>i and a constant fci >> 1 such that 

Vfc > ki, 



(3.26a) 
(3.26b) 



hm ej; = 



by ([3^ and ([3^. 

Now, it follows from (3.20), (3.26a), and the g-monotonicity of the numerical 
operator F (or G) that for k > max{/co, fci}, 

= F{6l ,^^Uh^{^k - hk),Sl^^^Uh,{^k),Sl^h^Uh^i^k + hk),^k) 
G{SlM^Uh^{£,k),dl^h^,UhdS.k),^k) 
> G{Sl^f,^Uh,i^k) + ek,Sl^f,^Uh^{^k),Ck)- 

Thus, by ( |3.18[ ), ( |3.26b| ), the consistency of P (or G), and ( |3.22p we get 

= liminf F (5^ ;,^ u/ij^ - (^^ji.Wh, (&), (^^ + /ife),Cfc) 

AC— >-00 

= liniinf G((5^,,^u/i,(Cfc) +efc,'5^^;,^Mhfc(6),6) 

> F*( lim Sl f^ Uh^{S,k),xo) 

fc->oo 

= (liminf (5^ i^u{xo), xq) 

> F*(limsup(5^ ,jU(a;o), a^o) 

where we have used the fact that F^, is decreasing in its first argument to obtain the 
last two inequalities. This is true by the definition of F* and Definition |3.1[ 

Step 3: We consider the general case ip G G^(f2) which is alluded in Step 2. Recall 
that u — if is assumed to have a local maximum at xq. Using Taylor's formula we 
write 

(p{x) = Lp{xo) + (p^{xq){x ~ xo) + ^(Pxx{xq){x - xq)^ + o{\x - XoP) 
p{x) + o{\x - xop). 
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For any e > 0, we define the following quadratic polynomial: 



ip{xo) + ipxixo){x - xo) + 



e + 



{x - xoY 



Trivially, p%^{x) = 2e + (px:x{xo) and ip{x) -p^ix) = o{\x - xa\'^) - e{x - xq)'^ < 0. 
Thus, if — p^ has a local maximum at xq, Therefore, u — p'^ has a local maxi- 
mum at xq. By the result of Step 2 we have Ft{p%^{xo) ^ xq) < 0, that is, F*(2e + 
fxxixo),xo) < 0. Taking liminfe_j.o and using the lower semicontinuity of Fa, we 
obtain > liminfe_j.o i^*(2e + (pxxixo), xq) > Fa,{Lpxxixo), xq). Thus, u is a viscosity 
subsolution of (3.1 )-( [33| . 

Step 4- By following almost the same lines as those of Step 2 and 3, we can 
show that if u — ip takes a local minimum at G f2 for some (p e C'^{fl), then 
F*{ipxx{xo)jXo) > 0. Hence, u is a viscosity supersolution of (3.1)-(3.3). 



Step 5: By the comparison principle (see Definition 3.3), we get u < u on fl. 
On the other hand, by their definitions, we have u < u on ft. Thus, u ~ u, which 
coincides with the unique continuous viscosity solution u of (3.1)-(3.3). The proof is 
complete. □ 

4. Two types of g-monotone finite difference methods. In this section we 
first construct two classes of practical finite difference methods of the form ( 3.4 ) . Using 
the first class of methods as examples, we then go through all the steps for verifying 
the assumptions of Theorem |3.4[ in particular, to present a fixed point argument for 
verifying the admissibility and stability. 

4.1. Finite difference methods with explicit numerical moments. We 

propose the following family of schemes with numerical operators: 

(4.1) Fi3{pi,p2,P3,x) := F{l3ipi + /32P2 + PzPz,x) + a{pi - 2p2 +P3), 

where {/3j}^^]^ are nonnegative constants satisfying /3i + /?2 + /^s = 1, and a is an 
underdetermined positive constant or function. 
Some specific examples from this family are 



(4.2) 
(4.3) 
(4.4) 



PI+P2+ P3 



,x] +a[pi - 2p2 +P3), 



Fiipi,P2,P3,x) := F 



F2ipi,P2,P3,x) := F{p2,x) +a(pi - 2p2 +P3), 



F3iPi,P2,P3,x) := F 



'Pi + 2p2 



P3 



a 



[pi - 2p2 + ps) ■ 



fact 



Remark 3. The term a(pi — 2p2 +P3) is called a numerical moment due to the 



26lU, + 5lU,+, 



a central difference approximation of u^xxxixj) scaled by . 
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4.2. Finite difference methods without explicit numerical moments. 

Given Pi,P2,P3 € R, let I{pi,P2iP3) denote the smallest interval that contains Pi,P2 
and p3, that is, 



I{Pi,P2,P3) ■■= [min{pi,p2,P3},niax{pi,p2,P3}]- 

Our first method in this family is the following Godunov type scheme (cf. 
the references therein). Its numerical operator F4 is defined by 



and 



(4.5) 
where 

(4.6) 



F4:{Pl,P2,P3,x) 



ext F{p,x), 

p6/(pi,P2,P3) 



ext 
pe/(pi,p2,P3) 



mm 

pG/(pi,P2,P3) 

max 

— . pe-f(pi,P2,P3) 



mm 

pi<p<P2 

min 

P3<p<P2 



if p2 > max{pi,p3}, 
if p2 < min{pi,p3}, 
if Pl <P2 < P3, 
if P3<P2 < Pl- 



Our second method in this family is a slight modification of the previous scheme, 
and its numerical operator, F5, is defined by 



(4.7) 
where 

(4.8) 



F5iPi,P2,P3,x) -.^ extr F{p,x), 



extr 

pe/(pi,P2,P3) 



mm 

pe-f(pi,P2,P3) 

max 

pe/(pi,P2,P3) 

max 

P2<P<P3 

max 

P2<p<pi 



if p2 > max{pi,p3}, 
if P2 < min{pi,p3}, 
if Pl <P2 < P3, 
if P3 <P2 < Pl- 



It is not hard to check that both F4 and are consistent and g-monotonc nu- 
merical operators. 

4.3. Verification of consistency, g-monotonicity, admissibility and sta- 
bility for scheme (4.1). In this subsection we use the methods with numerical 
operator Fp as examples to demonstrate all the steps for verifying the assumptions 
of the convergence theorem. Theorem 3.4 As mentioned before, the consistency and 
g-monotonicity are easy to verify, but the verification of the admissibility and stability 
are more involved. For simplicity, we only consider the case that F is differentiable 
and there exists a positive constant 7 > such that 



dF 

> -1/7 >7r^ 
op 



-1- 



(4.9) 
Recall that 

Fi3{pi,P2,P3,x) := F{(iipi + I32P2 + /^sPs,^^) + "(pi - 2p2 +P3), 
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where /3i, /32 and ^3 are nonnegative constants such that /3i + /32 + /?3 = 1- 

TriviaUy, Fj3(p,p,p^x) = F{p,x). Hence, i^^ is a consistent numerical operator 
for each set of /3i,/32 and (see Remark [l] (c)). To verify the g-monotonicity, we 
compute 

OFb ^ dF dFp ^ dF „ dF^ ^ dF 

Opi op Op2 op Ops op 

Then Fp is g-monotone if 

|^>o, ?fi<o. |§?>a 

dpi dp2 Op3 

On noting that ^ < 0, solving the above system of inequalities yields 

dF 

(4.10) a > -max{/3i,^3}^. 

Thus, we have proved the following theorem. 

Theorem 4.1. Fj^ is g-monotone provided that 

(4.11) a > max{/3i,/33}7 



for J defined by (4.9). 

Next, we verify the admissibility and stability of the schemes. To this end, we 
consider the mapping A4p : U ^ U defined by 

(4.12) 6lU, = SlU, + pFp{SlUj^i,SlU„6lU,+i,x,), j=2,3,--- ,J-l. 



Let U := iU2,U3-r-- ,C/j_i)^ and U := (C/2,C/3,--- ,Uj-i)^. Then ( |4l2| ) can be 
rewritten in vector form as 

(4.13) AIJ ^ AV + pG{V), 

where A stands for the tridiagonal matrix corresponding to the difference operator 
SlUj and G(U) = (G2(U,X2),G3(U,X3),- •• ,Gj_i(U,xj_i))^ with 

G,{lJ,x,) ^ Fp{dlU,.i,dlU„6lUj+i,x,), j = 2,3,--- ,J-1. 

Mp is said to be monotone if U is increasing in each component of U. 

Proposition 4.2. Suppose that Fp is g-monotone, that is, (4.11) holds. Then 
the mapping Ai p is monotone for sufficiently small p > 0. 

Proof. Consider the following system 

(4.14) W.^SlU,, j = 2,3,... ,J-1, 

(4.15) W, = W, + pFp{W,-i,W,,W,+i,x,), j = 2, 3, • • • , J - 1, 

(4.16) 6lUj^W„ j = 2,3,--- ,J-1. 

Let -.U^W, Mf^ -.W^W, and X^^) :W -^U. Then, it is easy to verify 

that Mp can be written as a composition operator of tM'^^ A^(^) and Ai^"^\ that is, 
Mp ~ M^^^ o Mf^ o M'-^l 
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Since A is positive definite, so is v4 ^. Thus, both M'--^^ and M^-^^ are monotone 
in the sense that they preserve the natural ordering of t°°{Th)- Moreover, since 



dWj _ dFp dWj _ dFp 



dFf, 



dWj+i ^ dp3 



then the g-monotonicity of Fp imphes that 



9W, 



> 0, 



— > 



and 



> 



provided that 
(4.17) 



0<p< [2a + /32/7] ^ 



Thus, Mf'' is monotone, so is Mp := TW^^^ o M^P o 



(4.17). The proof is complete. □ 



M.'^^\ provided that p satisfies 



Theorem 4.3. Under the assumptions of Proposition 4-2 the finite difference 
scheme (3.4) with F = Fp is admissible and stable. 

Proof. By the definition of G(U), we immediately have G(U + A) = G(U) for 
any constant A. Hence, Mp{lJ + A) = A^p(U) + A, and we have Mp commutes with 
the addition of constants. Together with the monotonicity of Mp, it follows that Mp 



is nonexpansive in £°°{Th) (see [8j). Hence (3.9) holds with C = max{|Ma|, |ub|}, and 
we have the scheme is stable. 

To prove admissibility of the scheme, let 



(4.18) SlV,^SlV,+pFp{SlV,^,,SlV,,SlV,+,,x,), j=2,3, 



, J-1. 



Subtracting (4.18) from ( |4.12[ ) and using the mean value theorem we get 



(4.19) 



1 + p 



dF, 



fl 

dp2 
ri 



Hence, 
(4.20) 



|U-V||,^ < (l + p[(/3i+/33)7-V7])l|U-V||,^ < -||U-V||,o., 



which holds for (/3i + /33) < 1/7^ and p > i [l /7 - (/3i + Ps)-/] . Thus, ( |420| imphes 



that the mapping Mp is contractive. By the fixed point theorem we conclude that 
Mp has a unique fixed point U, which in turn is the unique solution to the finite 
difference scheme (3.4) with F = Fp. The proof is complete. □ 

Remark 4. We note that the choice Pi — [3^ ~ ^ o-nd /32 = 1 trivially satisfies 
all of the restrictions in the proofs for any a > 0. We also note that the role of 
the numerical moment will be further explored numerically for degenerate elliptic test 
problems in section^ 
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5. Numerical Experiments. In this section, we perform a series of numerical 
tests to demonstrate the accuracy and the order of convergence for the various pro- 
posed numerical schemes. As before, we assume a uniform mesh. We use the Matlab 
built-in nonlinear solver fsolve for all tests, and, unless otherwise stated, we fix the 
initial guess [/'^'^^ as the linear interpolant of the boundary data. Also, all errors are 
measured in the L°° norm. 

For most tests, we record the results using Fi and F4. Unless otherwise stated, 
the results for all of the proposed Lax-Friedrichs-like operators are analogous and the 
results for all of the proposed Godunov-like operators are analogous, even though the 
analysis that prompts Remark [i] suggests F2 could be considered preferable to Fi and 
F3 . For most of the examples we observe quadratic rates of convergence to the viscosity 
solution for the Lax-Friedrichs-like schemes. For both classes of numerical operators 
we observe the lack of numerical artifacts that are known to plague the standard FD 
discretization for fully nonlinear problems. However, for the Godunov-like schemes, 
this phenomena typically presents itself through the fact that the nonlinear solver 
fsolve fails to find a root. Thus, while both classes of schemes support the selectivity 
of the discretizations, the resulting nonlinear algebraic system appears to be better 
suited for fsolve when using the Lax-Friedrichs-like operators. 

We begin with a simple power nonlinearity that has a C°° solution. 

Example 1: Consider the problem 

-ul^ + x^ = 0, -1 < a; < 1, 
u{-l) = -l/6, w(l) = l/6, 

with the exact solution u{x) = 

Using the linear interpolant of the boundary data as our initial guess and approx- 
imating u with the various schemes above, we obtain the computed results of Table 
|5.1| and Figure [5?T| 



h 


Fi ,a = 1.5 


F4 


L°° error 


order 


L°° error 


order 


l.OOOOe-01 


2.71e-02 




6.40e-02 




5.0000e-02 


5.10e-03 


2.41 


6.40e-02 


0.00 


2.5000e-02 


1.03e-03 


2.31 


6.40e-02 


0.00 


1.2500e-02 


2.33e-04 


2.14 


1.07C-03 


5.90 


6.2500e-03 


5.58e-05 


2.06 


2.12e-02 


-4.31 



Table 5.1 
Rates of convergence of Example 1. 



The schemes F2 and F3 exhibit similar behavior as Fi, and F5 exhibits similar 
behavior as F4. Thus, the Lax-Friedrichs-like schemes do exhibit a quadratic order 
of convergence as expected. On the other hand, the Godunov-like schemes converge 
inconsistently. This inconsistency is mostly due to fsolve failing to find a root. 

If we fix our initial guess as the approximation computed by Fi with a = 1.5 and 



h = 0.1, we get the results of Table 5.2 



Thus, Godunov-like schemes converge with high levels of accuracy when the non- 
linear solver has a sufficiently good initial guess. Since the Godunov-like schemes are 
very sensitive towards the initial guess for fsolve, it is hard to characterize a rate of con- 
vergence. We also observe in Table [5T] that the error for h = 0.1, 0.05, 0.025, 0.00625 
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EKample 1 witli scheme = 1 and N = 41 using alplia = 1 5 Example 1 with soheme = 4 and N = 81 




Fig. 5.1. Computed solutions of Example 1 



h 


Fi 1.5 


Fi 


L°° error 


order 


L°° error 


order 


l.OOOOe-01 


2.71e-02 




8.24e-08 




5.0000e-02 


5.10e-03 


2.41 


1.58e-06 


-4.26 


2.5000e-02 


1.03e-03 


2.31 


1.60e-05 


-3.34 


1.2500G-02 


2.33e-04 


2.14 


9.06e-05 


-2.51 


6.2500e-03 


5.58C-05 


2.06 


1.42e-02 


-7.29 



Table 5.2 
Rates of convergence of Example 1. 



is consistent with the error of the initial guess for the Godunov-hke schemes. In 
contrast, the Lax- Friedrichs- like schemes converge for a much wider range of initial 
guesses. 

The next example concerns the 1-D Monge- Ampere equation. 
Example 2: Consider the problem 

-M^^ -1-1 = 0, 0<x<l, 
u(0) = G, u(l) = 1/2. 

This problem has exactly two solutions 

JL- f \ I9 /\ -^2 

w (x) = -x", u (xj = — -x -I- a;, 

where is convex and u~ is concave. However, is the unique viscosity solution 
that preserves the ellipticity of the operator. 

Using U'^^'^ as the linear interpolant of the boundary data, the computed results 
with both types of schemes are given in Table |5.3| 

We note that the Lax-Friedrichs-like schemes converge to the unique ellipticity 
preserving solution (i.e., convex solution) for a > sufficiently large. However, if 
a < with \a\ sufficiently large, the Lax-Friedrichs-like schemes converge to u~ . The 
convergence to u~ for a < is expected since u~ is the unique solution that preserves 
the ellipticity of the PDE u^^ — 1 = 0. Forming the corresponding Lax-Friedrichs-like 
scheme and multiplying by —1 is equivalent to letting a < in the above formulation. 
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h 


Fi , a = 1 


Fi ,a = -l 


F4 


L°° error 


order 


L°° error 


order 


L°° error 


order 


l.OOOe-01 


2.54e-03 




2.54e-03 




1.17e-01 




5.000e-02 


6.36e-04 


2.00 


6.36e-04 


2.00 


1.21e-01 


-0.05 


2.500C-02 


1.59e-04 


2.00 


1.59e-04 


2.00 


1.24e-01 


-0.04 



Table 5.3 
Rates of convergence of Example 2. 




Example 2 with scheme = 1 and N = 21 using alpha = -1 



Example 2 with scheme = 4 and N = 21 




Fig. 5.2. Computed solutions of Example 2 with different parameter a. 



We also test the benefit of using a Lax-Friedrichs-like scheme as ojpposed to the 
standard 3-point finite difference method. We approximate u using F2 for varying 
values of a, using the linear interpolant of the boundary data as our initial guess. 
The computed results are given in Table [5^ 

We remark that letting a = corresponds to the standard 3-point finite differ- 
ence method, which does not converge in the above example. Instead, it behaves 
similarly to the Godunov-like schemes in that the nonlinear solver cannot determine 
a good direction to move from the initial guess. Thus, the Lax-Friedrichs-like schemes 
have a mechanism for giving the nonlinear solver a good direction towards finding a 
root. When a is sufficiently large, the schemes converge. When a is not sufficiently 
large, while the schemes may not converge, they have a tendency to move towards 
the correct solution. Furthermore, we can see that the Lax-Friedrichs-like schemes 
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h 


F2 , a = 6 


F2 ,a = 0.05 


F2 , a = 


L°° error 


order 


error 


order 


L°° error 


order 


l.OOOe-01 


3.07e-02 




1.18e-01 




9.00e-02 




5.000e-02 


8.51C-03 


1.85 


3.31e-02 


1.83 


1.15e-01 


-0.35 


2.500C-02 


2.14G-03 


1.99 


3.03e-02 


0.13 


1.15e-01 


-0.00 



Table 5.4 

Performance of a Lax-Friedrichs-like scheme with various a 



Example 2 with scheme = 2 and N = 41 using alpha = S Exannple 2 with scheme = 2 and N = 41 using alpha = 05 





D 0.2 0.4 O.B 8 1 



Fig. 5.3. Computed solutions by a Lax-Friedrichs-like scheme with various a 

converge quadratically for a bigger than the theoretical lower bound with only a 
small cost in the level of accuracy. Thus, when dealing with a problem that has an 
unknown optimal bound for a, large a values can be used. A shooting method for 
decreasing a allows the scheme to gain accuracy while maintaining the benefits of the 
Lax-Friedrichs-like schemes. 

If we first use Fi with a = 1 to approximate w on a coarse mesh with h = 0.1, 
and then we interpolate the result to get an initial guess for the two proposed schemes 
and the 3-point finite difference method, we get the results of Table |5.5[ Thus, we 
see that the Godunov-like schemes and the standard finite difference formulation now 
converge to with high levels of accuracy given a sufficiently good initial guess. In 
fact, they both converge to the same limit. 

To the contrary, if we use Fi with a = — 1 to approximate u on a coarse mesh with 
h = 0.1 and then interpolate the result as an initial guess, we obtain the results of 
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h 


Fi , a = 1 


F4 


F2 , a = 


L°° error 


order 




order 


L°° error 


order 


l.OOOe-01 


2.54e-03 




9.96e-15 




9.96e-15 




5.000e-02 


6.36e-04 


2.00 


4.54e-13 


-5.51 


4.54e-13 


-5.51 


2.500e-02 


1.59e-04 


2.00 


1.46e-10 


-8.33 


1.46e-10 


-8.33 


1.250e-02 


3.97e-05 


2.00 


9.85e-10 


-2.75 


9.85e-10 


-2.75 



Table 5. -5 

Performance of the standard 3-point scheme 



Example 2 with scfeme = 4 and N = 41 Example 2 with solieme = 2 and N = 41 using alpha = 




Fig. 5.4. Computed solutions by a Godunov-like scheme and the standard 3-point scheme 



Table 5.6 Clearly, none of the schemes converge to u^. Moreover, the Lax-Friedrichs- 
like schemes and the Godunov-like schemes do not converge to u~ even if is close to 
u~ . Instead, fsolve finds no solution when using the two proposed schemes. Thus, the 
Lax-Friedrichs-like schemes and the Godunov-like schemes appear to only consider m+ 
to be the solution of the PDE. Since u"*" is the unique viscosity solution of the PDE, 
lack of convergence to u~ for the Lax-Friedrichs-like schemes for a > sufficiently 
large and for the Godunov-like schemes is consistent with theory. 

In contrast, the standard 3-point finite difference method does converge to u~ . 
When given a sufficiently good guess, the 3-point finite difference method will converge 
to any one of the two solutions. Furthermore, the discretization can create artificial 
solutions that will attract the standard 3-point finite difference method. On the other 
hand, the monotonicity of our proposed schemes prevent the discretizations from 
having multiple solutions. 



h 


Fi , a = -1 


F4 


F2 , a = 


L°° error 


order 


J^OO gJ-J-QJ. 


order 


L°° error 


order 


l.OOOe-01 


2.68e-02 




2.56e-03 




2.24e-14 




5.000e-02 


5.61e-03 


2.25 


2.54e-03 


0.01 


8.82e-13 


-5.30 


2.500e-02 


1.26e-02 


-1.16 


2.54e-03 


0.00 


8.83e-12 


-3.32 


1.250e-02 


1.41e-02 


-0.17 


2.54e-03 


-0.00 


1.63e-09 


-7.53 



Table 5.6 

Performance of a Lax-Friedrichs-like scheme with a ■- 



The next two examples deals with Bellman type equations. 
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Example 2 wHh scheme = 1 and N = 1 1 using alpha : 




EKampIS 2 witli scheme = 2 and N = 1 1 using alpha = 




Fig. 5.5. Computed solutions of a Lax- Friedrichs- like scheme with a = — 1. 

Example 3: Consider the problem 



min {-AeUxx - S (x)] = 0, 



for 



-1 < a; < 1, 



if a; < 0, 
if a; > 0. 



This problem has the exact solution u{x) — x\x\ . We also note that this problem 
has a finite dimensional control parameter set. 

Using the linear interpolant as the initial guess, we obtain the results of Table 
|5.7[ We observe that the Godunov-like scheme converges and both schemes exhibit 
quadratic convergence for this example. 




h 




Fi 


L°° error 


order 


L°° error 


order 


l.OOOe-01 


1.29e-01 




9.60e-03 




5.000e-02 


4.67e-02 


1.46 


2.50e-03 


1.94 


2.500C-02 


1.46e-02 


1.68 


6.25e-04 


2.00 


1.250C-02 


4.18e-03 


1.80 


4.70e-01 


-9.55 


6.250e-03 


1.13e-03 


1.89 


4.72e-01 


-0.01 


3.125e-03 


2.95e-04 


1.93 


4.72C-01 


-0.00 



Table 5.7 
Rates of convergence of Example 3. 



Now we consider a Bellman problem with infinite dimensional control parameter 



set. 



Example 4: Let 

inf 

-1<8(2;)<1 



such that e _L°°([2,4]), and consider the problem 

2 < a; < 4, 



inf {-Bu^^ + 9'^u + x ^| = 0, 



u(2) = 4, u(4) = 16. 
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Example 3 with scheme = 1 and N = 321 using alpha = 1 Example 3 with scheme = 4 and N = 41 



Fig. 5.6. Computed solutions of Example 3. 



This problem has the exact solution u{x) = with the corresponding control 9{x) = 

Let the initial guess be given by the linear interpolant of the boundary data. 
Then, we obtain the results of Table |5.8[ 



h 


Fi , a = 0.5 




L°° error 


order 


L°° error 


order 


l.OOOe-01 


3.07e-01 




5.59e-01 




5.000e-02 


9.88e-02 


1.64 


4.96e-01 


0.17 


2.500e-02 


3.09e-02 


1.68 


5.10e+00 


-3.36 



Table h.i 
Rates of convergence of Example 



Both schemes have a hard time finding a root for h small, although the Lax- 
Friedrichs-like schemes do converge towards u for larger values of h. 
Now we choose the initial guess 



3 3 16 

14^ +y' 



a simple cubic function that satisfies the boundary conditions. Then, ||[/'^"^ ^u||/,oo([2_4]) 
0.94, and we get the results of Table [5^1 



h 


Fi ,a = 0.5 


Fi 


L°° error 


order 


L°° error 


order 


l.OOOe-01 


3.07e-01 




6.74e-10 




5.000C-02 


9.88e-02 


1.64 


7.04e-08 


-6.71 


2.500e-02 


3.09e-02 


1.68 


3.41e-09 


4.37 


1.250e-02 


9.02e-03 


1.78 


8.09e-08 


-4.57 


6.250e-03 


2.47e-03 


1.87 


9.44e-01 


-23.48 



Table 5.9 
Rates of convergence of Example . 
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EKample 4 with scheme = 1 and N = 1 1 using aipha =05 



u_guess 



E^ampie 4 with scher^e = 4 ahd N = 1 1 



Fig. 5.7. Computed solutions of Example 4- 



Thus, the Lax-Friedrichs-like schemes again converge with a rate of almost 2. 
Also, the Godunov-like schemes converge with high levels of accuracy for h > 0.0125, 
but for smaller h, fsolve fails to find a root . 

We remark that this problem can also be approximated by using a splitting al- 
gorithm. The operator can be split into an optimization problem for 9 and a linear 
PDE problem for u, and then a natural scheme is to successively approximate 9 and 
u starting with an initial guess for 9. For the above approximations, the nonlinearity 
due to the infimum was preserved inside the definition of the operator. 

The final example considers a problem whose solution is not classical. 
Example 5: Consider the problem 



-y-L + 8 sign(a;) = 0, 
'u(-l) = -l, u{l) = l, 



-1 < X < 1, 



with the exact solution u{x) — x\x\ G C^([— 1, 1]). 

Using the linear interpolant of the boundary data as the initial guess, we obtain 
the results of Table ISlOl 



h 


i^i , a = 1.5 


F4 


L°° error 


order 


L°° error 


order 


l.OOOe-01 


1.59e-02 




2.40e-01 




5.000e-02 


3.76e-03 


2.08 


2.50e-01 


-0.06 


2.500e-02 


9.40e-04 


2.00 


2.50e-01 


0.00 


1.250e-02 


2.35e-04 


2.00 


6.69e-06 


15.19 


6.250e-03 


5.88e-05 


2.00 


2.05e-01 


-14.90 



Table 5.10 
Rates of convergence of Example 5. 



We clearly see the quadratic rate of convergence for the Lax-Friedrichs-like schemes. 
The Godunov-like schemes only converge for h — 0.0125. For larger h, the scheme 
returns the initial guess after failing to find a root. For the test with smaller h, 
the scheme returns a slightly improved approximation after reaching the maximum 
number of iterations. 
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Fig. 5.8. Computed solutions of Example 4. 



If we fix our initial guess as the approximation formed by Fi with a — 1.5 and 



h — 0.1, we then get the results of Table 5.11 



h 


Fi , a = 1.5 


F4 


L°° error 


order 


L°° error 


order 


l.OOOe-01 


1.59e-02 




1.84e-08 




5.000e-02 


3.76e-03 


2.08 


4.05e-06 


-7.78 


2.500e-02 


9.40e-04 


2.00 


8.85e-06 


-1.13 


1.250e-02 


2.35e-04 


2.00 


6.50e-06 


0.45 


6.250e-03 


5.88e-05 


2.00 


7.78e-06 


-0.26 



Table 5.11 
Rates of convergence of Example 5. 



As observed in the previous examples, we see that the Godunov-like schemes 
converge quickly with high levels of accuracy, thus making it difficult to characterize 
a general rate of convergence. 

6. Conclusion. We have presented a new framework for constructing and ana- 
lyzing consistent, g-monotone, and stable finite difference methods. The newly pro- 
posed consistency and g-monotonicity criterion are not only simple to understand, 
but they are also easy to verify in practice. The key concept of the framework is the 
"numerical operator" , which plays the same role as the "numerical Hamiltonian" does 
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in the successful monotone finite difference framework for first order fully nonlinear 
Hamilton- Jacobi equations. To construct practically useful finite difference methods 
which can be easily implemented on computers, we have also presented a guideline 
for designing finite difference methods which fulfill the structure criterion of the pro- 
posed finite difference framework. The key concept in this regard is the "numerical 
moment" , which plays the same role as the "numerical viscosity" does in the successful 
finite difference framework for first order fully nonlinear Hamilton-Jacobi equations. 
Moreover, we gave some numerical evidences and argued that "numerical moments" 
provide an indispensable mechanism and ability for a finite difference scheme to be 
able to converge to the viscosity solution of the underlying second order fully nonlinear 
PDE problem. To a certain degree, the work of this paper bridges the gap between the 
state-of-the-art of finite difference methods for second order fully nonlinear PDEs and 
that for first order fully nonlinear Hamilton-Jacobi equations. Although the results 
of this paper are confined to the one spatial dimension case, they are also expected 
to hold in high spatial dimensions; that result will be presented in a forthcoming 
companion paper [15]. 
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